function [theta2,deScalingX,deScalingZ,hx_1,hx_2] = ...
    Var1Threshold_GMM_closedForm_macro(VarU_x2,Cov10U_x2,Cov01U_x2,...
    x2,macros,econAct,thresholdLevel,allowDeSaling,h0RegimeOn,hxRegimeOn)
% This function computes the closed form solution for the VAR(1) model 
% with threshold effects when the factors and macro display measurement errors

% Extending the output from th regression filter with macro variables
nm          = size(macros,1);
[nx2,T]     = size(x2);
factors     = [macros;x2];
nx          = size(factors,1);
VarU        = zeros(nx,nx,T);
Cov10U      = zeros(nx,nx,T);
Cov01U      = zeros(nx,nx,T);
for t=1:T
    VarU(nm+1:nm+nx2,nm+1:nm+nx2,t)   = VarU_x2(:,:,t);
    Cov10U(nm+1:nm+nx2,nm+1:nm+nx2,t) = Cov10U_x2(:,:,t);
    Cov01U(nm+1:nm+nx2,nm+1:nm+nx2,t) = Cov01U_x2(:,:,t);
end
% Defining regimes
regime_1    = econAct >= thresholdLevel; %a boom
regime_2    = econAct < thresholdLevel;  %a recression

% The process for the pricing factors
Xmat          = factors(:,2:T);
if h0RegimeOn == 1 && hxRegimeOn == 1
    a = [regime_1(1:T-1,1)';
        regime_2(1:T-1,1)';
        factors(:,1:T-1).*repmat(regime_1(1:T-1,1)',nx,1);
        factors(:,1:T-1).*repmat(regime_2(1:T-1,1)',nx,1)];
elseif h0RegimeOn == 0 && hxRegimeOn == 1
    a = [ones(T-1,1)';
        factors(:,1:T-1).*repmat(regime_1(1:T-1,1)',nx,1);
        factors(:,1:T-1).*repmat(regime_2(1:T-1,1)',nx,1)];
elseif h0RegimeOn == 1 && hxRegimeOn == 0
    a = [regime_1(1:T-1,1)';
        regime_2(1:T-1,1)';
        factors(:,1:T-1)];
elseif h0RegimeOn == 0 && hxRegimeOn == 0
    a = [ones(T-1,1)';
        factors(:,1:T-1)];
end             
VarUa         = zeros((1+hxRegimeOn)*nx+1+h0RegimeOn,(1+hxRegimeOn)*nx+1+h0RegimeOn);
Amat          = zeros(nx,(1+hxRegimeOn)*nx+1+h0RegimeOn);
if h0RegimeOn == 1 && hxRegimeOn == 1
    for t=1:T-1
        VarUa(2+1:2+nx,2+1:2+nx)           = VarUa(2+1:2+nx,2+1:2+nx) + VarU(:,:,t)*regime_1(t,1);
        VarUa(2+nx+1:2+2*nx,2+nx+1:2+2*nx) = VarUa(2+nx+1:2+2*nx,2+nx+1:2+2*nx) + VarU(:,:,t)*regime_2(t,1);
        Amat(:,2+1:2+nx)                   = Amat(:,2+1:2+nx) + Cov10U(:,:,t+1)*regime_1(t,1);
        Amat(:,2+nx+1:2+2*nx)              = Amat(:,2+nx+1:2+2*nx) + Cov10U(:,:,t+1)*regime_2(t,1);
    end
elseif h0RegimeOn == 0 && hxRegimeOn == 1
    for t=1:T-1
        VarUa(1+1:1+nx,1+1:1+nx)           = VarUa(1+1:1+nx,1+1:1+nx) + VarU(:,:,t)*regime_1(t,1);
        VarUa(1+nx+1:1+2*nx,1+nx+1:1+2*nx) = VarUa(1+nx+1:1+2*nx,1+nx+1:1+2*nx) + VarU(:,:,t)*regime_2(t,1);
        Amat(:,1+1:1+nx)                   = Amat(:,1+1:1+nx) + Cov10U(:,:,t+1)*regime_1(t,1);
        Amat(:,1+nx+1:1+2*nx)              = Amat(:,1+nx+1:1+2*nx) + Cov10U(:,:,t+1)*regime_2(t,1);
    end
elseif h0RegimeOn == 1 && hxRegimeOn == 0
    for t=1:T-1
        VarUa(2+1:2+nx,2+1:2+nx) = VarUa(2+1:2+nx,2+1:2+nx) + VarU(:,:,t);
        Amat(:,2+1:2+nx)         = Amat(:,2+1:2+nx) + Cov10U(:,:,t+1);
    end
elseif h0RegimeOn == 0 && hxRegimeOn == 0    
    for t=1:T-1
        VarUa(1+1:1+nx,1+1:1+nx) = VarUa(1+1:1+nx,1+1:1+nx) + VarU(:,:,t);
        Amat(:,1+1:1+nx)         = Amat(:,1+1:1+nx) + Cov10U(:,:,t+1);
    end
end
XXinv                = (a*a'-VarUa)\eye(size(VarUa,1));
ha                   = (Xmat*a'-Amat)*XXinv;
Wxhat                = (Xmat-ha*a)';
if h0RegimeOn == 1 && hxRegimeOn == 1
    h0_1                 = ha(1:nx,1);
    h0_2                 = ha(1:nx,2);
    hx_1                 = ha(1:nx,2+1:2+nx);
    hx_2                 = ha(1:nx,2+nx+1:2+2*nx);
elseif h0RegimeOn == 0 && hxRegimeOn == 1
    h0_1                 = ha(1:nx,1);
    h0_2                 = ha(1:nx,1);
    hx_1                 = ha(1:nx,1+1:1+nx);
    hx_2                 = ha(1:nx,1+nx+1:1+2*nx);    
elseif h0RegimeOn == 1 && hxRegimeOn == 0
    h0_1                 = ha(1:nx,1);
    h0_2                 = ha(1:nx,2);
    hx_1                 = ha(1:nx,2+1:2+nx);
    hx_2                 = ha(1:nx,2+1:2+nx);
elseif h0RegimeOn == 0 && hxRegimeOn == 0    
    h0_1                 = ha(1:nx,1);
    h0_2                 = ha(1:nx,1);
    hx_1                 = ha(1:nx,1+1:1+nx);
    hx_2                 = ha(1:nx,1+1:1+nx);
end    
VarWxHat             = zeros(nx,nx);
for t=1:T-1
    omegaHat = VarU(:,:,t+1)...
        + hx_1*VarU(:,:,t)*regime_1(t,1)*hx_1' ...
        + hx_2*VarU(:,:,t)*regime_2(t,1)*hx_2' ...
        - Cov10U(:,:,t+1)*regime_1(t,1)*hx_1' - hx_1*Cov01U(:,:,t+1)*regime_1(t,1) ...
        - Cov10U(:,:,t+1)*regime_2(t,1)*hx_2' - hx_2*Cov01U(:,:,t+1)*regime_2(t,1);    
    VarWxHat = VarWxHat + 1/(T-1-2*(nx+1))*Wxhat(t,:)'*Wxhat(t,:) - 1/(T-1)*omegaHat;
end

% The process for economic activity
z                            = econAct(2:T,1)';
b                            = [ones(1,T-1); econAct(1:T-1,1)'; factors(:,1:T-1)];
varUb                        = zeros(2+nx,2+nx);
varUb(2+1:2+nx,2+1:2+nx) = sum(VarU(:,:,1:T-1),3);
ZZinv                        = (b*b'-varUb)\eye(2+nx);
hb                           = (z*b')*ZZinv;
% Compare to ordinary regression!
%[B,BINT,R,RINT,STATS] = regress(z',[ones(T-1,1) econAct(1:T-1,1) factors(:,1:T-1)'])
%STATS(1,1) gives the R2
% [B,BINT,R,RINT,STATS] = regress(z',[ones(T-1,1) econAct(1:T-1,1) factors(3:4,1:T-1)']);
% hb = [B(1,1) B(2,1) 0 0 B(3,1) B(4,1)]
%[B,BINT,R,RINT,STATS] = regress(z',[ones(T-1,1) econAct(1:T-1,1) factors(2:4,1:T-1)']);
%[B,BINT,R,RINT,STATS] = regress(z',[ones(T-1,1) econAct(1:T-1,1) factors(1,1:T-1)'.*factors(3,1:T-1)'+factors(4,1:T-1)'.*factors(4,1:T-1)' ]);
%hb = [B(1,1) B(2,1) B(3,1) 0 0 0 0]
%WzHat = R;
%[B,BINT,R,RINT,STATS] = regress(z',[ones(T-1,1) econAct(1:T-1,1)])
%STATS(1,1) gives the R2

WzHat                        = (z-hb*b)';
gama0                        = hb(1,1);
gamaz                        = hb(1,2);
gamax                        = hb(1,2+1:2+nx)';
VarWzHat                     = 0;
for t=1:T-1
    VarWzHat = VarWzHat + 1/(T-1-(nx+2))*WzHat(t,1)^2 - 1/(T-1)*gamax'*VarU(:,:,t)*gamax;
end

deScalingX = 1;
[SwxHat,checkWx] = chol(VarWxHat,'lower');
[SwzHat,checkWz] = chol(VarWzHat,'lower');
if allowDeSaling == 1
    idx = 1;
    while checkWx ~= 0 && idx <= 100
        % Reducing deScaling
        deScalingX   = deScalingX*0.0;
        XXinv                = (a*a'-VarUa*deScalingX)\eye(size(VarUa,1));
        ha                   = (Xmat*a'-Amat*deScalingX)*XXinv;
        Wxhat                = (Xmat-ha*a)';
        if h0RegimeOn == 1 && hxRegimeOn == 1
            h0_1                 = ha(1:nx,1);
            h0_2                 = ha(1:nx,2);
            hx_1                 = ha(1:nx,2+1:2+nx);
            hx_2                 = ha(1:nx,2+nx+1:2+2*nx);
        elseif h0RegimeOn == 0 && hxRegimeOn == 1
            h0_1                 = ha(1:nx,1);
            h0_2                 = ha(1:nx,1);
            hx_1                 = ha(1:nx,1+1:1+nx);
            hx_2                 = ha(1:nx,1+nx+1:1+2*nx);
        elseif h0RegimeOn == 1 && hxRegimeOn == 0
            h0_1                 = ha(1:nx,1);
            h0_2                 = ha(1:nx,2);
            hx_1                 = ha(1:nx,2+1:2+nx);
            hx_2                 = ha(1:nx,2+1:2+nx);
        elseif h0RegimeOn == 0 && hxRegimeOn == 0
            h0_1                 = ha(1:nx,1);
            h0_2                 = ha(1:nx,1);
            hx_1                 = ha(1:nx,1+1:1+nx);
            hx_2                 = ha(1:nx,1+1:1+nx);
        end
        VarWxHat    = zeros(nx,nx);
        for t=1:T-1
            omegaHat = VarU(:,:,t+1)...
                + hx_1*VarU(:,:,t)*regime_1(t,1)*hx_1' ...
                + hx_2*VarU(:,:,t)*regime_2(t,1)*hx_2' ...
                - Cov10U(:,:,t+1)*regime_1(t,1)*hx_1' - hx_1*Cov01U(:,:,t+1)*regime_1(t,1) ...
                - Cov10U(:,:,t+1)*regime_2(t,1)*hx_2' - hx_2*Cov01U(:,:,t+1)*regime_2(t,1);
            VarWxHat = VarWxHat + Wxhat(t,:)'*Wxhat(t,:)/(T-1-2*(nx+1)) - deScalingX*omegaHat/(T-1);
        end
        [SwxHat,checkWx] = chol(VarWxHat,'lower');
        idx = idx + 1;
    end
    if checkWx ~= 0
        SwxHat = NaN(nx,nx);
    end
    
    idx = 1;
    deScalingZ = 1;
    while checkWz ~= 0 && idx <= 100
         % Reducing deScaling
        deScalingZ   = deScalingZ*0.0;
        ZZinv        = (b*b'-varUb*deScalingZ)\eye(size(VarUb,1));
        hb           = (z*b')*ZZinv;
        gama0        = hb(1,1);
        gamaz        = hb(1,2);
        gamax        = hb(1,2+1:2+nx)';
        VarWzHat     = 0;
        for t=1:T-1
            VarWzHat = VarWzHat + 1/(T-1-(nx+2))*WzHat(t,1)^2 - 1/(T-1)*gamax'*deScalingZ*VarU(:,:,t)*gamax;
        end
        [SwzHat,checkWz] = chol(VarzWHat,'lower');
    end
    if checkWz ~= 0
        SwzHat = NaN;
    end
end

% Saving out
theta2 = setTheta2_threshold_macro(h0_1,h0_2,hx_1,hx_2,SwxHat,gama0,gamaz,gamax,SwzHat);

end

